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Dynamo theory offers the most promising explanation of the generation of the suns 
magnetic cycle. Mean field electrodynamics has provided the platform for linear and 
nonlinear models of solar dynamos. However the nonlinearities included are 
(necessarily) arbitrarily imposed in these models. This paper conducts a systematic 
survey of the role of nonlinearities in the dynamo process, by considering the behaviour 
of dynamo waves in the nonlinear regime. It is demonstrated that only by considering 
realistic nonlinearities that are non-local in space and time can modulation of the basic 
dynamo wave be achieved. Moreover this modulation is greatest when there is a large 
separation of timescales provided by including a low magnetic Prandll number in the 
equation for the velocity perturbations. 
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1. INTRODUCTION 

Magnetic activity is found in many astrophysical bodies on all length- 
scales. It is now widely believed that stellar magnetic fields are 
generated by hydromagnetic dynamos in which the field is maintained 
against ohmic dissipation by the velocity fields within the star. 
Research into dynamo theory has received much attention in the 
recent past and many new approaches have recently been introduced. 
Specific attention has often focussed on the origin of the solar 
magnetic field, as the sun is the only star in which the magnetic field 
can be observed in detail. 
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The basic eleven year solar cycle is modulated on a longer timescale. 
For approximately 70 years during the seventeenth century, the 
amplitude of magnetic activity was dramatically reduced (Eddy, 1976; 
Ribes and Nesme-Ribes, 1993). This minimum in activity (The 
Maunder Minimum) can also be detected in the abundance records 
of terrestrial ,4 C and 1(J Be, with both records indicating the occurrence 
of several previous minima in magnetic activity (Stuiver, 1994; 
Beer et ai s 1991) 

Numerical studies of the generation of magnetic field in the sun have 
usually fallen into one of three categories (for more details see Weiss 
and Tobias, 1997). Numerical integration of the fully self-consistent 
dynamo equations (Gilman, 1983; Glatzmaier, 1985) demonstrate that 
dynamo cycles can be found. However these integrations are expensive 
to perform and are only feasible using over-simplified physics. 
Consequently, a number of alternative mean-field models have been 
studied (see e.g., Brandenburg et ai y 1989; Kitchatinov et a/., 1994) 
that attempt to give a greater insight into the temporal and spatial 
behaviour of nonlinear dynamo models. These models have met with 
varying degrees of success, producing modulated cycles and chaotic 
solutions, but the exact role of the various competing nonlinear 
mechanisms (e.g., quenching) is unclear. Moreover the modulation 
found in these models is different to the behaviour found in the sun. 
Modulation of the basic cycle is associated with large changes in the 
parity (symmetry) of solutions, whilst the sun undergoes significant 
change in the amplitude of oscillation with only very small changes in 
parity (Watari, 1996; Yau, 1988). 

The third approach is to use the recent advances made in the field of 
nonlinear dynamics to describe the behaviour of dynamo models. To 
illustrate how various aspects of solar magnetic activity may arise, it is 
often necessary to study toy models. These systems are not predictive, 
but they do enable the generic properties common to fully self- 
consistent dynamos and idealised mean field models to be easily 
understood (Weiss, 1993). This has been dubbed Astromathematics by 
Spiegel (1995). One approach is to concentrate on the alternation 
between active and quiescent periods. This behaviour, termed ‘on-off 
intermittency’ by Spiegel (1995), can be modelled in a number of ways. 
The first of these is aperiodic switching. This switching may rely on 
stochasticity [studied in detail in full PDE simulations by Choudhuri 
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(1992) and Ossendrijver and Hoyng (1996)] which is only likely to be 
of importance if the dynamo is marginal or, more naturally, 
deterministic chaos (Platt, 1994; Platt et al ., 1994; Spiegel, 1995). 
Indeed some models have been constructed that use both mechanisms 
(Schmitt et al ., 1996). Another technique is to consider the symmetries 
inherent in the problem and to construct a model that relies on the 
interactions (and resonances) between modes with different symme- 
tries. Knobloch and Landsberg (1996) find that for such a model two 
types of minima may be found. A complementary technique is to 
consider truncated representations of the mean-field dynamo equa- 
tions. For example, Weiss et al. (1984) extended Parker’s (1955, 1979) 
linear dynamo wave calculations to include nonlinear effects. They 
constructed a sixth order system of ordinary differential equations. 
This system, a complex generalisation of the Lorenz equations, yielded 
three successive Hopf bifurcations as the dynamo number D was 
increased. Multiply-periodic solutions lay on a three-torus in phase 
space. As D was increased further the torus was destroyed and 
aperiodic solutions were found. The model can be simplified to a fifth- 
order system whose bifurcation structure has been studied in detail 
(Jones et al 1985). These low-order models are open to the criticism 
that any behaviour found may be an artifact of the severe 
approximations made in truncating the system. For a model to be of 
astrophysical relevance, any properties should remain robust under 
small changes in the approximations and assumptions made. Covas 
et al. (1997) investigated the structural stability of a 'dynamic a-effect’ 
dynamo based on a model proposed by Schmalz and Stix (1991). They 
demonstrated that small changes in the level of truncation of dynamic 
forcing could have drastic consequences on the results of integrations. 
A different viewpoint was adopted by Tobias et al. (1995). They 
considered a robust third-order model constructed using the results of 
normal form theory and reproduced the bifurcation sequence found in 
the truncated model of Weiss et al. (1984). Although such results are 
expected to be robust, it has still to be demonstrated that this 
bifurcation structure would continue to be found in more realistic 
(PDE) models. 

To describe the dynamic behaviour of the dynamo it is vital first to 
understand the role of the nonlinear terms in the equations. In many 
mean-field models the nonlinear effect of the Lorentz force is included 
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to limit the growth of the magnetic instability. The back-reaction of 
the magnetic field on the velocity is often parametrised in such models. 
The small-scale motions, which lead to the a-effect and turbulent 
diffusivity, are sensitive to changes in the magnetic field strength 
(microdynamic quenching). An ad-hoc approach usually used in 
nonlinear mean-field models is to parametrise the action of the 
magnetic field on the trubulent motions by including a-quenching (see 
e.g., Stix, 1972). For consistency the turbulent diffusivity should be 
treated in a similar manner, because both tensors are derived from the 
same turbulence as first noted by Roberts and Soward (1975). The 
exact form of the quenching of the a-effect and turbulent diffusion is 
an open question requiring the investigation of sophisticated 
numerical models. What is clear however is that both of these effects 
should be suppressed by the action of the Lorentz force. 

The Lorentz force must also modify the large-scale velocity (a 
macrodynamic nonlinearity). This has often been included in models 
by quenching the differential rotation (w-quenching). More sophisti- 
cated models include the back-reaction of the small-scale magnetic 
field on the A-effect in the equation for the mean velocity field (e.g., 
Kitchatinov et a /., 1994). However the most natural consequence of 
the Lorentz force is the velocity perturbation driven by the mean 
magnetic field. This macrodynamic process has been termed the 
"Malkus-Proctor effect” after the work of Malkus and Proctor (1975), 
where it was used as the nonlinearity in an a 2 -dynamo model. Because 
the inclusion of this nonlinearity requires the solution of another 
equation coupled to the dynamo equation it has not been as 
extensively studied as a- or ^’-quenching. 

In this paper the effects of including both the microdynamic 
quenching mechanisms and the Malkus-Proctor effect on the proper- 
ties of dynamo waves are investigated. It is necessary to isolate the role 
of nonlinearities in the dynamo process from other effects that may be 
of importance in dynamo models. For example boundary conditions 
(Worledge et al . , 1997 Tobias et al ., 1997) and inhomogeneities 
(Meunier et al . , 1997) can both play a crucial role in determining the 
nature of solutions. The consequence of imposing various nonlinea- 
rities on wave-like solutions should therefore be determined before an 
attempt is made to include them in more realistic models. In this paper 
the simplest physically relevant model of a solar dynamo is considered. 
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The model is an extension of the one proposed by Parker (1993) for a 
dynamo acting at the base of the convection zone and characterised by 
a layer of strong shear lying beneath a layer where the a-effect is 
strong. In the Section 2 the important results of Parker are 
summarised and the features which are to be included in the new 
model are introduced. In Section 3 the new dynamo model is derived 
from the au> dynamo equations, and the linear theory for the new 
model is included in Section 4. Simple nonlinear behaviour is discussed 
in Section 4 where the consequences of including (microscopic) 
quenching mechanisms in the model are outlined. It is demonstrated 
that 17-quenching plays a very different role to simple equilibrating 
mechanisms such as a-quenching. The most important results are 
contained in Section 5 where the “Malkus-Proctor effect” is 
introduced and interesting temporal behaviour can be found. The 
bifurcation structure of the system is analysed and related to that 
found in earlier models. A discussion of the consequences of the results 
for solar and stellar dynamo theory is then included in Section 6. 


2. THE PARKER SURFACE DYNAMO WAVE 

The linear model considered by Parker (1993) demonstrates how a 
dynamo at the interface of the solar convection zone and radiative 
interior might work. The model provides the motivation this paper 
and, for this reason, a brief summary of the set-up and results for this 
system is included in this section. 

Parker constructed an extension of his initial (1955) ‘Dynamo Wave’ 
model in order to take into account the available observational data. 
A model of the base of the convection zone and the overshoot region 
which lies below is considered. Helioseismology (e.g., Thompson et ai, 
1996) indicates that there is a thin layer of radial shear at the base of 
the convection zone. This region of strong shear in the tachocline 
makes it a good candidate for the site of production of toroidal flux. 
However magnetic buoyancy of this strong azimuthal field needs to be 
suppressed in order to prevent the flux escaping to the surface on a 
timescale too short for the effective operation of a dynamo. How 
might the confinement of the toroidal field to this region be achieved? 
Parker argues that the total amount of azimuthal flux beneath active 
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latitudes is 10 23 - 10 24 Maxwells. If all this flux were to be concentrated 
into a thin region at the base of the convection zone the field intensity 
would need to be of the order of 2.5 kilogauss or more. The strong 
field would act back on the turbulence, suppressing the eddy diffusivity 
(via a nonlinear mechanism) and confining the field to a thin region. 
The problem is now that a field of this strength would also quench the 
a-effect, meaning that there would now be no generation of poloidal 
field in this region. The solution is to notice that, if the field is strongly 
confined, the suppression of the turbulent diffusivity and a-effect is 
local. The toroidal field in the bulk of the convection zone is weak and 
hence the cyclonic turbulence would only be weakly modified. Indeed 
Charbonneau and MacGregor (1996) have shown that this type of 
interface dynamo can work even if the a-effect is as strongly 
suppressed as has recently been suggested (Vainshtein and Cattaneo, 
1992; Cattaneo and Hughes, 1996). The convection zone is therefore 
the favoured site for production of poloidal flux from toroidal flux and 
the dynamo cycle is completed. The toroidal flux generated by the 
strong shear in the tachocline diffuses slowly into the convection zone 
where it is converted by the cyclonic turbulence into poloidal flux. This 
flux may then diffuse back to the overshoot region where it can act as a 
source for production of toroidal field. Here the region of generation 
of toroidal and poloidal fields are separated in space and the dynamo 
has to rely on diffusion to operate. For this reason it is less efficient 
than models which have the a-effect and w-effect in the same location. 
Other models with the generating regions in different places (e.g., 
Steenbeck and Krause, 1969; Kleorin and Ruzmaikin 1981) have 
previously been proposed. 

The model is constructed as follows. Parker considers local 
Cartesian coordinates with the origin at the base of the convection 
zone [as shown in Fig. 1(a)]. The system is two-dimensional and the 
surface z = 0 is the upper boundary of a region (z < 0) of uniform shear, 
q - ch'o/dz , where the eddy diffusivity is reduced to a value n and there 
is no a-effect. The region z>0 is convective and hence has a large 
turbulent diffusivity (ry) and a strong a-effect, but no shear; a diagram 
of the cartesian model is shown in Figure 1(b). Parker considers plane- 
wave solutions proportional to to the linear dynamo equations in 
the 2 layers, and the critical parameter is the non-dimensional dynamo 
number D = Ga/rfk 3 where k is the wavelength of the solution in the 
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CONVECTION ZONE 
STRONG TURBULENCE 


OVERSHOOT REGION 
STRONG SHEAR 


z =— 1 


FIGURE I The local Cartesian coordinates used by Parker (1993). 

x-direction. For D < D crit all plane wave solutions decay away as 
/— > oo. However for D>D criu R(<r)>0 and exponentially growing 
solutions may be found. Parker examines two cases in detail. The first 
is when the ratio of the diffusivities in the two layers is unity so 
(j 2 = njr] = 1 . In this case the dispersion relation that determines a 
reduces to a simple form and Z) cr i t may easily be found. This simple 
calculation establishes the existence and nature of the solutions. They 
take the form of a dynamo surface wave propagating along the 
interface between the region of cyclonic convection and the region of 
shear. The second case considered is the limit n 77 so that /r <§C 1. 
This represents the situation of very small flux transport in the region 
of differential rotation (z<0). Here the calculation yields a vertical 
scale for the toroidal field in z < 0 that is substantially smaller than that 
in z>0 and the azimuthal field in z<0 is confined to a thin layer 
pressed up against the underside of z = 0. The field in z>0 is more 
diffuse. Further calculation yields the interesting result that in order 
for the calculation to be consistent, the relationship 

nb 2 m . M <XT]B 2 mM ( 1 ) 

must hold where 6 max and £ max are the maximum values of the 
toroidal field in the regions z<0 and z>0 respectively. The relation- 
ship in (1) is kinematic, and for a given /1 2 determines the ratio 
This scenario is therefore entirely self-consistent if the 
eddy-diffusivity is suppressed inversely with the magnetic energy 
density. 
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Because the calculation is linear it does not contain any information 
about how the magnetic field may act in the nonlinear regime. Parker 
hypothesised that some form of dynamical balance is achieved in the 
sun, and used the example of the action of Maxwell stresses on Alfven 
waves to illustrate a possibility for this balance. In this paper we 
examine how the addition of nonlinearities changes the surface wave 
model, and which properties of the linear model remain robust in the 
nonlinear regime. In the next section, the nonlinear dynamo problem 
is discussed and the results of previous nonlinear models are presented. 


3. DERIVATION OF THE MODEL 


The starting point of the derivation is the mean field dynamo 
equation 

(9B 

— = Vx(ux B + QB- 77 VX B), (2) 


where B and u are the mean magnetic field and velocity field of the 
plasma respectively. The turbulent diffusivity (77) and the a-effect arise 
from the mean interaction of the small-scale fluctuating magnetic and 
velocity fields. The dynamo equations are considered in the two- 
dimensional Cartesian coordinates considered by Parker. Immediate 
progress can be made by decomposing the magnetic field into toroidal 
and poloidal parts, i.e. 

B = B(x,z) y + V x [A(x,z)y], (3) 


Equation (2) is clearly linear in the magnetic field B and for a given 
velocity field solutions are either exponentially growing or decaying. 
As noted in the introduction, nonlinearities may be introduced in the 
model by parameterizing the back-reaction of the magnetic field on the 
turbulence (a and 77 -quenching). 

In order to include the effect of the magnetic field in driving a 
velocity perturbations (macrodynamic nonlinearities) it is necessary to 
consider the equation for the large-scale velocity u. The momentum 
equation for the mean velocity field u is given by 



+ (u- V)u 


= —Vp + V - <r + 2H x u + j xB + (4) 
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where p is the density, p is the pressure ST is the rotation of the star, <r is 
the stress tensor and j= l//x 0 V x B is the large-scale current. The 
effects of rotation and magnetic field on the velocity are clear in (4), 
and the other forces (e.g., gravity, thermal forcing, Reynolds stresses 
from the small-scale turbulence) are contained in the tensor . 

Progress is made by assuming that the large-scale velocity consists 
of two separate components. The first component would be present in 
the absence of a magnetic field, and is generated by purely 
hydrodynamic effects. This velocity is driven by rotational, gravita- 
tional and thermal effects and provides the basic state of the model. 
The second component arises owing to the back-reaction of the 
magnetic field via the Lorentz force. It is important to stress here that 
the Lorentz force generated by the dynamo will drive a meridional 
flow in addition to the toroidal velocity considered here. Here we 
ignore this effect and concentrate solely on the toroidal velocity 
perturbations. It will be interesting in the future to study models where 
the latitudinal shear is also included and to investigate the interaction 
of the driven meridional flow with this shear. 

In the model, the above statements correspond to setting 

u = U + v = (J{x, z) y + v(x, z) y, (5) 

where U is the imposed basic large-scale velocity [given by 
dVjtlz = G(z);G(z) is the shear defined by the Parker model] and v 
is the velocity perturbation that is driven by the macroscopic magnetic 
stresses due to the Malkus-Proctor effect. 

In the calculations that follow, the vertical extent of the layer is 
considered to be small compared to a density scale height and 
therefore the density p is set to be a constant. This is a reasonable 
assumption for a dynamo located in a thin region at the base of the 
solar convection zone, where the density scale height is large. Here we 
consider solutions of the coupled system (2) -(5). These equations are 
non-dimensionalised by setting 

x = Lx, / = L 2 /r)oi , r] = r}oh(z) , a = ao/(z), G(z) - ^\)g(z), 

(5a) 


v(z) - Lu; 0 v(z), A{z) = LB 0 <*qA(x,z), B(x,z) = B 0 B(x,z), (5b) 
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where L is a characteristic lengthscale in the z-direction. In the limit 
where the shear is strong in comparison with the a-effect (the au; limit) 
the following equations are obtained 


dA_ d 1 A d 2 A 
~dt^ h d^ + ~dJ 


+fB , 


(6a) 


dB 

~dl 


\d 2 B 

L 

d 2 B 

w 

dz\ 


dhdB „ 

+ + D 

dz dz 


g{z) + 


dv 


dA dv dA 

dx dx dz T 


(6b) 


9V h 


d 2 v d 2 v 
dx 2 + d ? 


dh dv 

+ T aiS + T 


dBdA 

dz dx 


dA dB 
dz dx 


(6c) 


where r = uo/r]o(the magnetic Prandtl Number) is assumed to be a 
constant and 


D = aou^ r= _^aoL 
Vo pofioVo^o 

Here Z), the non-dimensional dynamo number, is a measure of 
the efficiency of the dynamo giving the strengths of the o and u; 
effects. The parameter T is therefore the product of the ‘a-effect 
Reynolds number’ R a = aoL/rjo and the (modified) Elsasser Number 
A = B\j pofitoVowo- It is important to note at this point that, because 
this parameter contains the quotient ao/tuo, the sign of T is the same 
as that of the dynamo number Z), i.e. if the dynamo number is chosen 
to be negative because the product R a R^ is negative then for consist- 
ency the parameter T is negative also. 

Equation (6) is to be solved in a two-dimensional layer -1 < z < 1, 
0 < x < -V m ax for suitable choices of the functions j\ g, and h. Here we set 

f(z, B)= If 5[ , + tanh ^( z - z 3)][l“ e_25( '” 2o) ^ ]Qo{B), z > z 0 , 
f(z)Qo{B) } ~ \ Q 

( 8 ) 


g(z,B) =\ 
g(z)Q K (B)j 


5 [1 -tanh £’(z-z 2 )](l -e 25{: Z|), ]e K (fl), z<n, 

0, z>Z], 

(9) 
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h(z,B) -rfo - ) _ ( i[l+tanh£(z-z 3 )][l-c 25(z 2o)2 ] 0 M ( 5 ), z > z 0 , 

A(z)0#i(*) i i 7/o, 

( 10 ) 

where £*, z 0 , 2^22,23 are all real parameters that may in principle be 
varied. The nonlinear quenching mechanisms to be employed are given 
by QaiB ), Q k {B) and Q^{B), with the level of quenching depending on 
the value of <7, k and /1 respectively. Both local and global (where the 
quenching is inversely proportional to the average energy over the 
period of a wave) quenching laws will be investigated. 

It is difficult to visualise the form of the a-effect, shear and turbulent 
diffusion from the mathematical formulae (8) -(10), and so an example 
of the imposed profiles for the a-effect, w-effect and turbulent 
diffusivity for a particular choice of parameters is shown in Figure 2 . 
This figure shows that the profiles are similar to those considered by 
Parker in his linear model (shown in Fig. 1 ). The a-effect is limited to 
the convection zone (2 > 0) and the shear is strongest in the overshoot 



FIGURE 2 Profiles for the shear (dashed), alpha-effect (solid) and diffusivity (dot- 
dashed) showing the z-dependence. Here zo = -0.75, z\ — 0.0 and Z 3 = 0.75 are fixed. 
The parameters chosen are E = 10, Z 2 — 0.5. Notice the region of overlap of the a-effect 
and shear. 
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region (z<0). The dependence of solutions on the details of these 
profiles is discussed below. 


4. LINEAR THEORY 

4.1. Derivation of Equations and Boundary Conditions 

The aim of the linear calculation is to determine the critical dynamo 
number, i.e., the value of D at which instability to dynamo action sets 
in. The perturbed velocity v only arises as a result of nonlinear 
interactions, and so for the linear magnetic theory v = 0. Moreover the 
nonlinear quenching is not considered and so Q a = Q K = = 1. In 

light of the work of Parker (1955,1993) the initial bifurcation is 
expected to be a Hopf bifurcation to travelling wave solutions; the 
preferred direction of travel depends on the sign of the dynamo 
number. Travelling wave solutions of the form A = a(z)e ikx+(Tt , 
B = b(z)e tkx+(7t are sought where a = A + iuj. Thus the dynamo waves 
have frequency uj and have a growth-rate A. the problem now reduces 
to finding the value of D = D t . for which A = 0. We set a(z) = a R (z) + 
iai(z),b(z) = bfi(z) + ibj(z), substitute into (6) and equate real and 
imaginary parts to obtain the set of coupled ODEs 

h{z)b" R (z) + h l {z)h f R {z) - [/i(z)/r -f A]/^(z) — -ub/ + kDg(z)a /, 
h{z)b;{z) T ti(z)h'j(z) - [h(z)k 2 + \]b/(z) - wb R - kDg(z)a R , 

h{z)a” R {z) - \h(z)k 2 + A }a R (z) = -ua, - f(z)b R , ' 1 

h(:)a,{:) - [h{z)k 2 + A )a,(z) = ua R -f(z)h,. 

Initially the boundary conditions 

B/4 

A = B = 0 at z= -l, — =fi = 0atz=l (12) 

are chosen as reasonable conditions at the bottom and top of the layer. 
Physically these conditions correspond to setting B v = B v = 0 at the top 
boundary and B = 0 at z=-l. It is desirable to show that changing 
these boundary conditions has only a weak effect on the main 
properties of the solution. The linear solutions are periodic in x with 
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period The linear theory will determine the preferred scale (i.e., 
wavenumber) in the .v-direction and this will be used in the nonlinear 
regime to select the aspect ratio of the box. The boundary conditions 
imply that 


b R = b, = a R = ci[ = 0, z = - 1, 
b R = bi = a'# =■ c/j — 0, z = -j- 1 . 

The equations and boundary conditions constitute a two-point 
boundary value eigenvalue problem that may be solved numerically. 
Two numerical schemes were employed, the first used the program 
NRK (see Gough et at., 1976) which employs the Newton-Raphson- 
Kantorovitch method using a second order centred-difference dis- 
cretisation. The results were checked using the NAG routine 
D02HBF, which uses a shooting method to solve the eigenvalue 
problem. 

4.2. Properties of Eigensolutions 

The profiles shown in Figure 2 are selected as a realistic scenario for 
the base of the convection zone. They are a slight modification of the 
model proposed by Parker (1993). It seems likely that the transition 
between the region dominated by shear and that in which there is an 
efficient a-effect will be fairly sharp, but that this will be a smooth 
transition (rather than a discontinuous jump). Moreover, it would be 
surprising if there were no overlap of these two regions, as the 
convection may overshoot into the region of shear (and indeed the 
shear may persist slightly into the base of the convection zone). It will 
be shown later that the important properties of the solutions are 
unaffected by the precise choice of profile, and it is only the form of the 
eigenfunctions that change significantly. 

There is the symmetry D— ► — D, x — > —x inherent in the problem, 
and for this reason only the case D < 0 need be considered. The results 
for D >0 are the same but with the direction of travel of the dynamo 
wave reversed. To determine the preferred scale in the x-direction, the 
critical dynamo number for a given horizontal wavenumber k is 
calculated. The calculation is then repeated for other values of k and 
the curve of marginal stability is shown in Figure 3. The curve has a 
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a) b) 




FIGURE 3 (a) Graph of critical dynamo number, - D c , against wavenumber, k , for 
the parameter values chosen in Figure 2. THe curve has a minimum and therefore there 
exists a preferred lengthscale for the problem, (b) Graph of frequency u against 
wavenumber, k for the same parameter values as in (a) (with D = D C ). The frequency is a 
monotonically increasing function of wave number and tends to a limit as k — > 0. 


minimum, and therefore a preferred scale. The minimum dynamo 
number is \D ( \ = —D* — 40.24 and this occurs at k = k* =0.78. This 
value suggests that it is sensible to consider the region 0 < x < 2n in the 
nonlinear regime which corresponds to a wavenumber k= 1 for which 
the critical dynamo number is Z) r = -42.0. In a similar way the second 
eigenvalue, corresponding to the frequency oj of the travelling wave at 
the onset of dynamo action, can also be calculated as a function of 
wavenumber, and this curve is shown in Figure 3(b). Here a; is a 
monotonic function of wavenumber tending to infinity as k is 
increased and to a finite limit, a; ->1.46 as k— >0. For k= 1 the 
frequency is u; = 2.50, and so the period of the wave is given by T= 
2tt/M = 2.51. 

A contour plot of one wavelength of the eigenfunction with k= 1 is 
shown in Figure 4(a, c). The toroidal field is generated by the shear in 
the region z < 0.5 and diffuses to the region of a-effect (z > 0) where the 
poloidal field is regenerated. For this reason the bulk of the toroidal 
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a) b) 





FIGURE 4 (a) One wavelength of the eigenfunction for the toroidal field (B) in the 

linear problem, with parameters chosen as in Figure 2. (b) Graph of (B 2 ) versus z, which 
clearly shows the z-dependence of the linear eigenfunction. The toroidal field is generated 
by the u/-eflect from the a-effect and is strongest in, and just below, the region of overlap 
of the a-effect and u;-effect. (c) As (a) but showing contours of the vector potential A. (d) 
As (b) but for (A 2 ) versus z. The graph shows that {A 2 )j{B 2 ) is small (although only the 
relative size of the fields is determined by the linear theory). 


field is generated and stored in (and just below) the region of overlap 
of a-effect and shear. Similarly the poloidal field generated by the 
a-effect from the toroidal field lies in z> 0, in (and just above) 
the overlap region as in Figure 4(b, d) where the z-dependence of the 
eigen functions is shown. These results are consistent with those found 
by Parker (1993) for his linear model. 

The calculations were then repeated for values of the parameters 
that change the shape of the profiles for the a-effect, shear and 
diffusivity, given in equations (8) -(10). These results are summarised 
in Table I. As expected, the larger the region of overlap of the shear 
and a-effect the more efficient the dynamo; therefore the critical 
dynamo number is smaller for profiles that have large regions of 
overlap. It is also apparent that the preferred wavenumber does not 
vary significantly when the parameters are changed. All the calcula- 
tions undertaken have k= 1 as the preferred integer wavenumber. This 
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TABLE 1 This table shows the minimum critical dynamo numbers, -£>*, and the 
corresponding wavenumbers, A + , for various parameter choices for the linear problem 
defined in (1 1) and (13). Also shown are the values of the critical dynamo numbers and 
frequencies for k = 1 These will be the critical dynamo numbers in the nonlinear 
problem. It is clear, that although changing the profiles has a qualitative effect on these 
values, the critical wavenumbers and frequencies are all similar. The critical values for 
E 10, z 2 = 0.5 are shown in bold type. These are the parameters that will be used in the 
nonlinear theory 


E 

2l 

Effect of Changing n, u> and t\ profiles 

D* k* —D at k = I 

uj at k = I 

1 

0.5 

9.59 

0.58 

12.0 

1.07 

10 

0.5 

40 24 

0 78 

42.0 

2.50 

20 

0.5 

55.19 

0.83 

56.7 

3.02 

50 

0.5 

67.33 

0.87 

69.1 

3.65 

1 

0 

12.97 

0.58 

16.3 

1.03 

10 

0 

110.0 

0.75 

116.9 

1.41 


would seem to indicate that the chosen scale for the eigenfunction 
arises not from any lengthscale that may be introduced by the profiles 
(i.e., the length of the overlap region or a typical length for the switch- 
off of the a-effect), but from the nature of the Parker model and the 
radial extent of the domain. 

The frequency is a measure of the efficiency of dynamo action. In 
the linear regime a more effective measure of the efficiency is the 
growth-rate of the solutions (and consequently the critical dynamo 
number), but in the nonlinear regime where solutions grow until 
equilibration a useful measure of efficiency is how quickly the dynamo 
can generate field of the opposite sign (i.e., how quickly the field can 
reverse) and this is given by the frequency of oscillation. As noted by 
Parker (1993), a dynamo with regions of a-effect and shear separated 
in space is less efficient than one where the poloidal and toroidal fields 
are generated in the same location, because of the time taken for 
the toroidal field produced by the shear to diffuse to the region where 
the poloidal field is generated (and vice-versa). Larger separation of 
the two regions implies less efficient dynamo action, and it should 
therefore also be the case that the dynamo is made less efficient if the 
diffusivity between the two regions is reduced - the magnetic field 
then takes longer to diffuse between them. This explains the changes in 
the frequency shown in Table I - the frequency is higher (indicating a 
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more efficient dynamo) for profiles that have a larger region of overlap 
and a larger turbulent diffusivity in that region. This explanation will 
be useful in analysing the nonlinear results when the back-reaction of 
the magnetic field on the profiles (and hence the frequency) is 
considered. 

The effects of changing the boundary conditions of the linear 
problem and introducing anisotropies into the turbulent diffusion on 
the nature of the solutions are now investigated. Solutions are found 
for the same equations with modified boundary conditions on A at the 
top of the box, so that (11) is now solved subject to the boundary 
conditions 


bR = bj = cir — ai = 0, z — — 1, 

(14) 

b R = bj — a R = a f = 0, z = -f 1 . 

The parameter values defining the profiles for the shear, a-effect and 
turbulent diffusivity were chosen to be the same as in the first 
calculation. It is found that there are no qualitative differences in the 
curves of marginal stability- the only differences are quantitative. For 
this calculation D* = -74.05 and this occurs at k = k* = 1.09. Clearly 
k= 1 will again be the preferred mode in the nonlinear calculation. 
(The slight increase in the size of the critical dynamo number can be 
attributed to the fact that the .4 = 0 boundary condition at z= 1 makes 
poloidal fields harder to excite near the boundary.) Also the frequency 
is as before a monotonic function of wavenumber. 


4.3. Summary of Linear Results 

The dependence of critical dynamo numbers and wavespeed on 
wavenumber has been calculated for a variety of parameter values and 
boundary conditions. The results clearly show that the preferred 
lengthscale is not determined by the precise choice of parameters in the 
model and that other properties of the solution are robust and are not 
sensitive to the choice of either profile or boundary conditions. Rather 
they are determined by the general characteristics of the model (large 
turbulent diffusivity, shear and a-effect largely separated). Therefore, 
these properties should also apply in the nonlinear regime and indeed 
in more elaborate systems based on the same model. 
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5. QUENCHING MECHANISMS 

In this section the consequences of including simple quenching 
mechanisms in the model are investigated and a systematic survey of 
these types of nonlinearities is undertaken. The most usual quenching 
mechanism employed in dynamo models is a-quenching where the 
magnetic field acts back on the small-scale helical motions. It has been 
argued (e.g., Vainshtein and Cattaneo 1992; Cattaneo and Hughes 
1996) that this effect can lead to suppression of the growth of magnetic 
field at levels much lower than equipartition. In this section no 
assumptions are made a priori about the strength of this effect -rather 
it is considered to be a parameter of the model. Similarly u>quenching 
and j ]- quenching (representing modification of the angular momen- 
tum profile and turbulent diffusion by the magnetic field) are 
considered. 


5.1. Non-Local Quenching 

The simplest extension of the model to the non-linear regime is to 
consider simple quenching mechanisms acting in isolation. In the past 
these have been the preferred non-linearities in many dynamo models, 
and so the results should act as a benchmark for the more complicated 
non-linearities considered later. Here the Malkus-Proctor effect is 
ignored by setting T = 0 in (6c). Now (6c) decouples and (6a, 6b) may 
be solved in isolation with v set to be zero. 

The most basic forms of quenching to consider are those where the 
a-effect, ^-effect and turbulent diffusion are quenched inversely with 
the average magnetic energy of the wave. That is we set 


Qo = 


l 


1 +cr £ 2 




1 


1 + kB 2 


e„ = 


\+nB 2 


(15) 


where again an overbar indicates a horizontal average. As discussed 
extensively in Tobias (1996b), the nature of the solution inherent in the 
choice of periodic boundary conditions makes this form of quenching 
a natural approximation. Local (in both x and z) quenching 
mechanisms will give qualitatively similar results as discussed in 
Section (5.2). 
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Now the a and uj effects and turbulent diffusion are not only 
functions of position but are also modified by the average toroidal 
field B 2 . Nonlinear travelling wave solutions of the form 

m 

B = y~]frj(z)siny(x - ct ) + bj(z)cosj{x - ct), 

7=1 (16) 

m v ' 

A = y^a)(z)siny(-v - ct) + a C j(z)cosj{x - ct), 

7=1 

can be found using a continuation method. Here c is the speed of 
propagation of the wave in the .v-direction and the wavenumber is set 
to unity. The Galerkin expansion is now substituted into (6a, 6b) to 
yield the system comprising the 4 m coupled nonlinear ordinary 
differential equations for the coefficients; i.e., for j = 1, m 

h[b” c (z) -fb c j(z)} + cjb](z) + Dgja](z) + h'bj(z) 
h[bf(z) -j 2 bj(z)\ - cjbj (z) - Dgja c j(z) + h'b'cj{z ) 
h[a" c {z) —j 2 dj(z)] + cja]{z) +/6'(z) 
h[a'J s (z) -j 2 a)(z)} - cja c j{z) +fb](z ) 

where ' denotes differentiation with respect to z. At this point it is 
important to recall that h , g, and / are all functions of both z and B 2 
and hence there is coupling between all the modes. This is apparent 
when it is noted that Parseval’s Theorem implies 



fi2(z) = £[&;(z) 2 +^(z) 2 ], 


= 2 ^(z)b l f(z) + b^z)bf(z)), 

(18) 


and so all the modes in the sum appear in every equation in (17). 

Equation (17) defines a new two-point boundary value problem, 
with the same boundary conditions as those initially considered in the 
linear theory (solutions are implicitly periodic at .v = 0, 2n). In this way 
both stable and unstable branches of travelling wave solutions to the 
problem may be located. 
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The stability of these solution branches is then determined by 
timestepping the equations using a pseudo-spectral discretisation in 
space together with an Adams- Bashforth timestepping scheme. The 
timestepping routine itself will only find stable solution branches but 
will also allow solutions other than travelling waves to be found. 

The nonlinear equations introduce three new parameters to the 
problem. These are a, k, and ji that represent the efficiency of the a-, 
u;- and //- quenching in the model respectively. The increase in the 
number of parameters in the problem means that an extensive search 
through parameter space is not feasible. However, the linear theory 
discussed earlier indicated that the results would not be too sensitive to 
particular selections in parameters that modify the profiles for the 
shear, o-effect and difTusivity. It therefore seems more sensible to 
explore the effect of the quenching parameters in the nonlinear regime, 
and to keep the ‘profile parameters’ fixed at £=10, z 2 ^ 0 (so the 
profiles are as in Fig. 2). The boundary conditions are also held fixed. 

5 A A. a- and (n-quenching 

The simplest extension to the nonlinear regime for the travelling wave 
solutions is to introduce a-quenching as an equilibration mechanism 
for the field (i.e., set a= 1 without loss of generality // = « = 0). The 
branch of travelling waves is followed from the initial Hopf 
bifurcation at D = — 42.0, as the magnitude of the order parameter D 
is increased. The branch is shown in Figure 5 where the generated 
magnetic field is a monotonically increasing function of the dynamo 
number D\ this is a feature of all the nonlinear solutions found in this 
section. The travelling waves remain stable for all the values of D 
investigated, which is also a result that holds for all further choices of 
parameters. The lack of subsequent bifurcations to more complicated 
time-dependent behaviour means that the effects of varying the 
parameters (and hence the nonlinearities) on the spatial structure and 
the wavespeed can be investigated with relative ease in this system. 
Clearly a different type of nonlinearity (other than simple quenching) 
is needed to produce more complicated time-dependence in this 
system. 

In Figure 5(b) the wavespeed c is plotted as a function of dynamo 
number. It can be seen that the wavespeed decreases slightly as —D 
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a) b) 




FIGURE 5 Q-quenching - nonlinear solution: (a) Average magnetic energy {B 2 ) at 
z = 0 versus dynamo number, —D for the case ct — 1, k = fi ~ 0. (b) Wavespeed c versus 
dynamo number with parameters as in (a). The wavespeed stays nearly constant 
(although a slight drop is noticeable). 


increases, so the waves slow down as \D\ is increased; at first sight this 
seems a strange result. The quenching mechanism acts so as to reduce 
the effective value of the dynamo number until D e rr ~ A The 
wavespeed should therefore remain approximately constant (as noted 
in Noyes et al 1984). The slight decrease in wavespeed is due to the 
modification of the profile for the a-effect by the local action ot the 
magnetic field, and can be understood in the context of the results of 
the linear theory discussed earlier. Figure 6 shows how the growth of 
magnetic energy in the area of overlap of the regions of a-effect and 
shear quenches the a-effect in that region. The a-effect is now confined 
to a thin region near the top of the layer, and there is no longer a large 
region of overlap. Hence the dynamo is less efficient (as discussed 
earlier) and the wavespeed is less than that for the critical value of the 
dynamo number. The nonlinear toroidal and poloidal fields have a 
similar spatial dependence as the eigenfunction of the linear problem, 
as shown in Figure (6a). This indicates that the quenching acts simply 
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(a) (b) 




FIGURE 6 Nonlinear solutions (a-quenching) — z-dependence: (a) Graph of ( B 2 ) 
(horizontal) versus z (vertical) for D= -1164, other parameters as in Figure 5. The 
toroidal field is concentrated at and just below z = 0 as in the linear eigenfunction, (b) 
Graph of o-effect (horizontal) versus z (vertical). The a-effect has been quenched to zero 
by the magnetic field throughout the bulk of the region, only reaching its unquenched 
level at z = 1 where the toroidal field is set to 0. 


as an equilibration mechanism for the field, and does not change the 
spatial dependence of the solutions. 

The case where k = 1 and // = a = 0 is now considered, so that only cu- 
quenching is operating. Figure 7 shows that again the magnetic energy 
at z = 0 is a monotonic function of D and the travelling wave solutions 
remain stable. Moreover the wavespread of solutions does not vary 
significantly as D is increased. In this case only a slight increase in 
wavespeed is observed for large dynamo numbers. The dynamo is 
therefore slightly more efficient for larger forcing. The solution for 
large dynamo numbers has a similar z-dependence to the eigenfunction 
of the linear theory, as in the case where a-quenching is considered in 
isolation. 


5.7.2. q-quenching 

So far the nonlinear quenching has acted solely as an equilibration 
mechanism but it is now necessary to include the effect of the Lorentz 
force on the diffusion properties of the turbulence. Now we set /z= 1, 
so that / 7 -quenching is in operation, (a-quenching also remains in 
operation as cr = 1 ) and it is expected that the z-dependence of the 
solutions will change as the turbulent diffusion is modified by the 
action of the magnetic field. 
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a) b) 




FIGURE 7 w-quenching nonlinear solutions: (a) As for Figure 5(a), but parameters 
are now k = 1, <7 = /z = 0. (b) As for Figure 5(b), but parameters are as in (a). Again the 
wavespeed c is largely unaffected (~25% change), although a slight increase in c may be 
detected. 


In Figure 8(a) the dependence of the magnetic energy on dynamo 
number is displayed and when this is compared with Figure 5(a) it is 
clear that the toroidal field strength (and hence magnetic energy) that 
may be created is significantly larger for the case when //-quenching is 
included. This effect can be explained using either mathematical or 
physical arguments. It is necessary to understand how the reduction of 
the local value of the turbulent diffusivity can lead to an increase in the 
strength of the magnetic field that may be generated by dynamo 
action. The mathematical argument is simply that the reduction of 
diffusion locally by //-quenching leads to an increase in the value of the 
effective local dynamo number D ]oc , This follows immediately from 
the definition 


D 


loc 


D 


h 2 (z,B 2 ) 


(19) 
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FIGURE 8 77 -quenching (together with a and tu-quenching): (a) As for Figure 5(a), but 
for parameter values <7 = k - 1. The magnetic energy is much larger for a given 
dynamo number than for the case with just q and ^/-quenching. (b) As for Figure 5(b), 
but with parameters in (a). THe wavespeed is reduced to nearly zero as the dynamo 
number is increased. This is a new effect introduced by the rj-quenching. 


and assuming that a- and u-quenching are kept at the same level. An 
increase in D ioc means that the local generation of field is more efficient 
for a given global dynamo number and this leads naturally to stronger 
fields. Physically, local dissipation via the turbulent convection is 
reduced by the action of the Lorentz force and this leads to the build 
up of magnetic energy locally and to significantly larger field strengths. 
At first there may appear to be a contradiction here, reducing 
dissipation leads to the build-up of stronger fields yet earlier it was 
argued that magnetic dissipation is vital to the Parker dynamo. How 
may the two be reconciled? The answer is that as the turbulent 
diffusion is quenched the form of the solution changes (see below) so 
that the regions of poloidal and toroidal field overlap. Diffusion is now 
not so essential to the operation of the dynamo as in the linear model 
because both the toroidal and poloidal fields are in the same location. 

The most striking effect of introducing 77-quenching can be seen in 
Figure 8 (b) where the wavespeed is plotted as a function of dynamo 
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number. In contrast to the cases where only a - and u>quenching were 
considered the wavespeed does depend significantly on the dynamo 
number. The wavespeed is a monotonically decreasing function of 
dynamo number, showing a sharp drop initially and appearing to tend 
to zero as D— >oo. This drastic reduction in wavespeed is due to the fact 
that the turbulent diffusivity has been almost completely quenched by 
the magnetic field and molecular diffusion is the only mechanism for 
the transport of magnetic flux. This makes the dynamo significantly 
less efficient as diffusion of poloidal field from the convection zone to 
the overshoot region and of toroidal field in the opposite direction are 
both suppressed. As the dynamo is less efficient the wavespeed c is 
decreased and the cycle period is increased for higher dynamo 
numbers. This result raises the question of how much the diffusivity 
may be quenched in astrophysical bodies -if the quenching is too 
strong the dynamo may cease to be efficient. 

Another indication that the quenching is no longer acting simply as 
an equilibration mechanism comes from the spatial dependence of the 
solution as shown in Figure 9. The z-dependence of the 77 -quenched 
solution (particularly of the toroidal field) is very different from the 
eigensolution. This is because there is no longer a sharp change in 
the diffusivity in the region where the magnetic field is generated, so 
the rapid change in preferred lengthscales that occurs at z = 0 in the 
linear theory (and in the nonlinear a- and u;-quenching models) is no 
longer present. As a result the field is no longer confined near to z = 0 
and can penetrate further in to the region z < 0 , where the shear is 
strongest and so toroidal field may be easily produced. The slanted 
appearance of the toroidal field is due to the action of the shear and 
the time taken for the toroidal field produced by the shear to diffuse 
upwards into the region z > 0. All these results demonstrate that the 77 - 
quenching mechanism has a very different effect on the nature of the 
solutions from the other types of quenching. It changes the form of the 
solutions and has a drastic effect on the wavespeed. 

5.2. Local Quenching 

The form of the (averaged) quenching selected in (15) is a simpli- 
fication consistent with the hypothesis that the solutions are local and 
therefore represent the magnetic behaviour at a fixed latitude. The 
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a) b) 



FIGURE 9 //-quenching-contour plots: (a) One wavelength of the toroidal field B< for 
the parameters in Figure 8 and £> = -1150. The toroidal field may now penetrate into the 
region r<0. (b) One wavelength of ,4, for the parameters as above. 


simplification to a one-dimensional (z-dependent) form for the 
quenching enabled the nonlinear solutions to be followed using a 
continuation method and the radial dependence to be studied in 
detail. However it is important to check that the main conclusions 
remain valid if the form of the quenching is changes to a iocaf (fully 
two-dimensional) form. Hence the calculations are repeated for the 
case where the quenching is local and the definitions in (15) are 
replaced by 

Qa= 1 + <tB 2 (.x, :) ’ Qk = 1 + kB 2 (x, 7 ) ’ Q " % + /,S 2 (.v, -)' 

( 20 ) 

The equation for the velocity perturbations is again decoupled by 
setting T = 0. The equations are now fully two-dimensional and are 
therefore too complicated to solve using a continuation method. 
Stable solutions are found using the time-stepping code described 
earlier. 
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The results of the integrations show that the solutions behave in an 
analogous manner to those for the case with averaged quenching. For 
example for simple local a-quenching only, we set cr= 1, k = /z = 0. For 
this choice of parameter values the travelling waves remain stable as D 
is increased into the highly nonlinear regime. As in the case where the 
quenching is proportional to the average magnetic energy of the wave, 
the magnetic field strength is a monotonic function of dynamo number 
and the speed of travel of the waves varies only marginally. For all 
other calculations where the velocity equation is decoupled and local 
quenching mechanisms are the sole nonlinearity, travelling waves 
remained stable. Hence the introduction of local dependence to the 
quenching models has very little effect on the nature of the solutions. 
As in the case of global quenching the temporal behaviour of solutions 
remains uncomplicated. 

Quenching mechanisms are particularly simplified parametrisations 
of the effect of the Lorentz force, being local in space and 
instantaneous in time. This inherent simplicity means that the form 
of the solutions is not sensitive to the exact choice of the nature of the 
mechanism i.e. whether the mechanism is averaged (one-dimen- 
sional) or local (two-dimensional). In Section 6 it will be demonstrated 
that this is not the case for the more-complicated (and realistic) 
Malkus-Proctor nonlinearity, where the particular form of the 
interaction can have large consequences for the nature of the solutions. 

5.3. Summary of Quenching Results 

In the above sections a thorough investigation of the nonlinear 
behaviour of the travelling waves has been conducted. The important 
results are summarised below. The travelling waves remain stable for 
all choices of parameters and all values of the dynamo number 
investigated. This indicates that the quenching mechanisms (whether 
local or averaged) considered do not naturally lead to modulation of 
the waves. It may be that a more realistic formulation of the problem 
where the Lorentz force is allowed to act back directly on the velocity 
will lead to more complicated time-dependence and this is investigated 
in Section 6. 

a- and ^-quenching have similar properties. Both act purely as 
equilibration mechanisms. They have little effect on the spatial 
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structure of the solutions. Also the speed of the travelling waves 
remains roughly constant as the dynamo number is increased. 
However, the introduction of 77 -quenching changes both the spatial 
structure of the waves and also the wavespeed. The wavespeed is 
reduced as the dynamo number (and hence magnetic energy of 
solutions) is increased. The dynamo becomes less efficient as the 
turbulent diffusion is reduced by the strong fields. The initial Hopf 
bifurcation to travelling waves may be subcritical if the ? 7 -quenching is 
allowed to be significantly stronger than the ^-quenching. 


6. THE MALKUS-PROCOTOR NONLINEARITY -ROUTES 
TO CHAOS 

In this section I will demonstrate that it is only by considering a 
sufficiently realistic (and therefore complicated) nonlinearity that a 
model can be constructed that exhibits transitions to complicated 
temporal behaviour. It is therefore vital to include the role of the 
Lorentz force in generating velocity perturbation via the Malkus- 
Proctor effect. These results will prove to be the most interesting and 
important when describing the modulation of the basic solar cycle. 

We therefore consider the fully-coupled nonlinear system given by 
(6). The boundary conditions that are imposed on the magnetically 
driven velocity v must be consistent with the (implied) conditions 
satisfied by the shear g. In this case, it is logical to choose periodic 
boundary conditions in the ^-direction and 

dv 

v(x,-l) = 0, -(*, 1) = 0, (21) 

as the conditions on the top and bottom of the box. 


6.1. Global (Averaged) Interaction 

6.1.1. Derivation of Equations 

The travelling wave models discussed so far are local models in the 
sense that they contain no latitudinal dependence of the parameters. In 



PROPERTIES OF NONLINEAR DYNAMO WAVES 


315 


this case it is also reasonable to adopt nonlinearities that are averaged 
over the period of a wave in the x-direction, as the primary interest is 
in the radial (z) dependence of solutions and their temporal behaviour. 
Calculations that include local quenching mechanisms have been 
considered, and local interaction via the Malkus-Proctor effect will be 
considered later in this section. For now, the a-effect, shear and 
turbulent diffusion remain quenched with the average magnetic field 
energy so that/(z, B 2 ), g(z, B 2 ) and h(z, B 2 ) remain as defined in 
(8) — (10). If the quenching is averaged in the latitudinal direction then 
for consistency the perturbation to the velocity should only arise due 
to the Lorentz Force averaged over the period of the travelling wave in 
the x-direction. In this case (6) becomes 


dA 

~di 

dB 

dt 


= h 


d 4 = rh 

dt 


d 2 A &A 
dx 2 ' dz 2 
d 2 B d 2 B 
dx 2 ~ dz 2 
d 2 v d 2 v 
dx 2 z 2 


+fl l 

dhdB n 

+ 7T7r + Z) 

dz dz 


/ \ dv 
8{Z)+ Jz 


dA 
dx ’ 


( 22 ) 


dhdv . . . 

+ r— ^- + sign (D) 
oz dz 


'dBOA 
dz dx 


dAdB 
dz dx 


where an overbar indicates a horizontal average over the period of a 
wave in the x-direction. These equations with the averaged effects of 
the Lorentz Force form a starting point for this investigation. 
Travelling wave solutions can be found by expanding the poloidal 
and toroidal fields in a Fourier expansion moving with the wavespeed 
as in (16). In a similar way ordinary differential equations for the 
coefficients aj, aj, b} and b] can be derived, though this time they are 
coupled to an equation for v(z). The equations become 

h[b'; c {z) -j 2 bj{z)] + cjbj{z) + D^g + J)X(Z) + h'bj(z) = 0, 

h[b'J s (z) -fb s j(z)\ - cjbj(z) - o(g + ^^(z) + h'b'f = 0, (23 ) 

h[af{z) - fa c j{z )] + cja*(z) +fb‘j(z) - 0, 
h[a"/{z) -j 2 a)(z)\ - cja^z) +fb](z) - 0, 
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for 7=1, m , together with 
d 2 \ d\> I 

Th d? + Th> Jz + Sign (Z)) 2 YsWi ~ b 'i a i - h ?j + h j“'fl = 0. (24) 

“ 7=1 

These equations (23) -(24), together with the appropriate boundary 
conditions on the toroidal and poloidal field defined in (12) and 


K- 1) = ^ (1) =°, 


(25) 


constitute a two point boundary value problem. The stability of these 
travelling wave solutions may be determined using the pseudospectral/ 
finite difference timestepping scheme described earlier. Indeed, as 
bifurcations to more elaborate time-dependent behaviour are ex- 
pected. this timestepping scheme will be useful in locating stable 
quasiperiodic and chaotic solutions. Here the properties of these 
travelling waves are investigated in the expectation that they will lead 
to some understanding of the Malkus-Proctor effect that may prove 
useful in analyzing the full problem. 


6.1.2. Properties of Travelling Waves: Wavespeed and Shear Profiles 

Rather surprisingly (and perhaps disappointingly), the travelling 
waves remain stable for all the combinations of parameter values 
selected. To obtain the required modulation therefore, one must 
consider the full problem, with local interaction of the magnetic and 
the velocity fields, defined by (6) (as in Section 6.2). The travelling 
waves do however show different behaviour than those found for 
simple quenching mechanisms. Initially the Malkus-Proctor effect is 
set to be the only nonlinearity in the problem and k, o and p are all set 
to be zero and the magnetic Prandtl number r is set to be unity. All 
other parameters relating to the shape of the profiles for the a-eifect, 
shear and turbulent viscosity are set to be the same as in Figure 2. 

Here (as before) negative dynamo numbers are studied. Figure 10 
demonstrates that both the magnetic energy and wavespeed are 
increasing functions of the magnitude of the dynamo number. 
The wavespeed increases by a significant amount. The z-dependence 
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FIGURE 10 Malkus-Proctor only - nonlinear solutions: (a) Average magnetic energy 
{ B 2 ) at z = 0 versus dynamo number, —D for the a — k = p = 0,r = 1 (b) Wavespeed c 
versus dynamo number with parameters as in (a). The wavespeed increases as the 
magnitude of the D is increased. 


of the toroidal field in the nonlinear regime is similar to that for the 
eigenfunction with the only change being that the Lorentz Force 
generates a shear tower down (i.e. at smaller values of z) than the 
region where the magnetic energy is concentrated. This shows one 
important difference between quenching mechanisms and the Malkus- 
Proctor effect -whilst quenching mechanisms may only act where the 
field is strong, the Lorentz force may drive shears in places with 
weak(er) fields. This is because the Lorentz Force term in the velocity 
equation depends on gradients of the field (and not soley on the 
absolute magnitude). If the other nonlinear effects are kept at zero, 
then changing the value of r in the averaged velocity equation has no 
effect on the nature of the solutions. Indeed it may be removed from 
the equations by rescaling the magnetic field. Hence if a = k — // = 0 
then t— 1 without loss of generality. 
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Having understood the case with the Malkus-Proctor effect only, we 
now wish to include the quenching mechanisms to see what effect (if 
any) they may have. A large number of combinations of parameter 
values are available to us, and here only a small subset is presented. 
The first quenching mechanism to be added is o-quenching The 
magnetic energy is less if both nonlinearities are included together than 
if the nonlinearities are considered in isolation (unsurprisingly). The 
dependence of the wavespeed is slightly more complicated; initially if 
both mechanisms are included the wavespeed is higher as the Malkus- 
Proctor mechanism tends to increase the wavespeed. However if D is 
increased lurther the speed drops below that of the case where the a- 
quenching is acting in isolation which demonstrates that the two 
mechanisms must interact in a more complicated manner than simply 
adding the two effects together. 

If a;-quenching and //-quenching are added the story is much the 
same - the magnetic energy is less if the Malkus-Proctor effect is 
included. For the wavespeed, the //-quenching is again the dominant 
nonlinearity, slowing the wave down for the reasons discussed in 
earlier sections - the Malkus-Proctor effect does have some effect in 
counteracting this behaviour, but it is clear that the efficiency of the 
dynamo is diminished by the suppression of the turbulent diffusion. 

6.2. Local Interactions 

The global interaction described in the previous secion failed to lead to 
the complicated time-dependence anticipated. For this reason we 
return to the full problem where the interactions of magnetic fields 
with velocity fields are local, and no averaging takes place. If the 
Lorentz force drives a large-scale velocity locally then for self- 
consistency the quenching mechanisms that parametrise the effects of 
the Lorentz Force on the small-scale turbulence should also depend on 
the local values of the magnetic fields strength. Hence the definitions of 
/, g and h representing the a-effect, shear and turbulent diffusion are as 
in the local models described in equation (20). 

The local quenching makes it very difficult to use continuation 
methods on the problem as the problem is now fully two-dimensional. 
For this reason the timestepping routine is used to find stable solutions 
for this problem. 
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6.2. L Malkus-Proctor Effect Alone 

As in the last section, the quenching mechanisms are initially switched 
off by setting <t = k = fi = 0 and the Prandtl number r is set to be unity. 
For \D\ < /) cri t = 42 both the magnetic field and velocity field decay 
away as predicted by linear theory but as the dynamo number is 
increased past this value dynamo action sets in at a supercritical Hopf 
bifurcation to travelling waves. For D just greater than this critical 
value stable travelling waves are found - the solution is periodic, but 
is constant in a frame of reference travelling with the wavespeed. 
Analysis of the solutions is easier if this basic oscillation associated 
with the speed of travel of the wave is removed, which can be 
accomplished by examining some global properties of the solutions. 
The average magnetic energy E B and the energy kinetic energy of the 
magnetically driven velocity field E K are calculated, where 

E b = ^ B 2 dxdz, E k = -t v 2 dxdz (26) 

where V is the domain of integration (0 < x < 2n,— 1 < 2 < 1). 

For travelling wave solutions we except E B — > const, and E K — ► const, 
and the trajectory to tend to a fixed point in the (E B , E K ) phase plane. 
For a weakly nonlinear solution both E B and E K are reasonably small. 

The solutions for the toroidal and poloidal field are shown in 
Figure 1 1 - these are very similar to the eigenfunctions, which is not 
surprising as the induced velocity v is very small and so will have only 
a weak effect on magnetic field generation. The figure shows the wave 
travelling along the interface from left to right; nearly a complete cycle 
of the wave, which has a period 2.5, is displayed. The wavespeed 
can therefore be calculated as 2.5 (1 d.p.) which is consistent with 
the linear theory. 

The magnetically-driven velocity perturbations are displayed in 
Figure 12. As the velocity is driven by the quadratic Lorentz force its 
x-dependence is different to that of the magnetic fields. In this case it is 
clear that B~B 0 sin (x — c7) and so it is natural that the velocity 
v~v 0 + V 2 sin 2{x—ct). It is the resultant shear that is important in 
equilibrating the dynamo, and it can be inferred that this is 
concentrated at small values of 2 . (The presence of this shear can be 
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FIGURE II Magnetic Field Travelling Waves. The figure shows the form of the 
solutions for A (top) and B (bottom) as time progresses for D = -50. The sequence 
should be viewed from left to right and top to bottom. Each snapshot is separated by a 
time t ^ 0.4992. Strong positive fields are shown in red whilst strong negative fields are in 
blue. (See Color Plate I). 


detected by the rapid change of colours from blue to red to green near 
—0.75.) In this respect this model is similar to the averaged model 
considered in the previous section as the driven shear appears at small 
values of z away from the region of maximum magentic field. This 
shows once again that the Malkus-Proctor effect is non-local and 
therefore is more sophisticated than simple quenching mechanisms. 
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FIGURE 12 Velocity Perturbations Travelling Waves. The figure shows the form of 
the solutions for r as time progresses for D —50. The velocity perturbations oscillate 
with half the period of the magnetic field. The snapshots correspond to those in 
Figure 1 1. The presence of a radial shear near z -0.75 can be inferred from the rapid 
colour change. (See Color Plate II). 


The velocity perturbations migrate in the same direction as the 
magnetic field with the same wavespeed, but because of the .v- 
dependence of the perturbations, the period of the travelling waves is 
half that of the magnetic fields; i.e. B => me 2/Lj '. 

If \D\ is increased further, the travelling waves grow in amplitude 
until they become unstable in a secondary Hopf bifurcation to 
modulated waves. For these modulated waves both the magnetic 
energy E H and the kinetic energy E K oscillate about a mean value as 
shown in Figure 13 with the period of oscillation being 0.46. Hence for 
these parameter values, trajectories in energy space tend to a limit 
cycle as shown in Figure 13(c). This limit cycle corresponds to a two- 
tor ns when the basic oscillation due to the travelling waves is 
reintroduced. 

The modulation is apparent in Figure 14 where the structure of the 
solution is now very different to that of the eigenfunction. After the 
secondary Hopf bifurcation the solutions not only travel but 'breathe' 
as they move along the interface. As well as this change in the 
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a) 




FIGURE 13 Modulated waves periodic energies, (a) Magnetic Energy /:# for D = 
-150. (b) Kinetic Energy E v for D~ -150. 


temporal behaviour of the solutions, the spatial dependence has clearly 
altered dramatically. The velocity field produced by the large-scale 
Lorentz force near z = — 1 is now strong enough to generate some 
toroidal field. Two properties of this induced field are immediately 
apparent from Figure 14. The first is that this field varies on a smaller 
scale than the field in r>-0.5 - in this layer Z?~sin 2x. This can be 
attributed to the fact that this field is generated by the velocity 
perturbation and not the imposed shear and that it is entirely 
contained in a region where the diffusion is small and so the natural 
lengthscale is also smaller. The second property concerns the speed of 
propagation of the waves. The magnetic field in this layer travels very 
slowly to the right - a fact that can be seen by comparing the position 
of the strong field in the first and last frame in Figure 14. The speed of 
propagation of this field can be estimated to be ^2^/27 = 0.23. The 
dynamo waves travel more slowly in this region because the diffusivity 
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FIGURE 14 Magnetic Field Modulated Waves. As for Figure 11, but £>=-150. 
The velocity field has acted back producing toroidal field at small values of z. This 
magnetic field has a smaller scale as it is produced by the velocity perturbations (which 
have a characteristic length-scale half that of the main field) in a region with a smaller 
value of the diffusivity. It travels slowly to the right for the same reason. (Sec Color Plate 
III). 

is small in this layer and the dynamo is therefore less efficient and 
hence slower. The magnetic field that is generated by the imposed 
shear for larger values of z continues to travel relatively quickly and is 
largely unaffected by the magnetic field in the layer of reduced 
diffusivity. There is some interaction however, when a fast-moving 
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region of strong positive field passes over a slow-moving region of 
positive field (as in the first snapshot on the left-hand side) these fields 
merge together and the radial (and latitudinal) gradients in the field 
will be small in the region. However, when strong positive field passes 
over strong negative field (as in the last snapshot on the right) large 
gradients in field are created. These are important in generating the 
velocity perturbations via the Malkus-Proctor effect. 

The velocity perturbations corresponding to this magnetic field are 
shown in Figure 15. The snapshots show that the velocity at the top of 
the domain resembles that for the travelling wave case. At the bottom, 
however, the velocity has quite a different form: four nodes, which 
oscillate and travel slowly to the right, appear in this layer of reduced 
diffusion. The size and stregth of these nodes at any given time can be 
understood in terms of the relative positions of the fast-moving and 
slow-moving magnetic field. If, as explained earlier, the proximity of 
strong negative toroidal field and strong positive field leads to large 
gradients in the field, then these gradients are responsible for the 
generation of velocity perturbations. 



HGURE 15 Velocity Perturbations Modulated Waves. The figure shows the form 
of the solutions for r as time progresses for D = -150. Velocity perturbations are drived 
by gradients in the magnetic field. (See Color Plate IV). 
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However this effect is not instantaneous - the velocity takes some 
time to respond to the magnetic field, and so there is a time-lag 
between the magnetic field and the generation of shear. It is well 
known from the work of Yoshimura (1975) that the introduction of a 
time-lag to the nonlinearities in the dynamo equations can lead easily 
to complicated time-dependence. The Malkus-Proctor effect, by 
including the equation for the large-scale velocity, introduces this 
time-lag in an entirely natural way. This is in contrast to the simple 
quenching mechanisms discussed previously, where the small-scale 
turbulent velocity field has an instantaneous response to the magnetic 
field which is local in space. 

If the dynamo number is further increased, the modulated waves 
persist, with an increase in amplitude and frequency. However for 
D = -330 the nature of the solutions is different as shown in Figure 1 6 
where the basic oscillation in the magnetic and kinetic energy is itself 
modulated on a longer timescale. This indicates that the solution has 
undergone a tertiary Hopf bifurcation where a third frequency of 



FIGURE 16 Quasipcriodic waves (a) As for Figure 13 (a), but Z)=-330. (b) As for 
Figure 13 (b), but D= —330. The energies are now quasipcriodic. 
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oscillation has been introduced to the problem. Motion now lies on a 
two-torus in energy-space (corresponding to a three-torus in the 
extended system where the basic wavespeed is included). The spatial- 
dependence for this type of solution is very similar to that for the 
modulated wave described earlier. 

As D is increased further this torus breaks down and the solutions 
become chaotic as shown in Figure 17. The ‘ghost' of a torus still 
remains for this parameter value (D= -400) as shown in Figure 17(c). 
In this case however the nature of the spatial dependence of solutions 
is somewhat different as shown in Figure 18. The solution is still 
characterised by two layers of interacting dynamo waves travelling at 
different speeds. The horizontal scale for the variation of the field is 
similar in both layers and the two layers interact substantially. This 
implies that the velocity perturbations are now' as important as the 
imposed shear at the top of the box. If D is increased further, chaotic 
solutions of ever-increasing complexity are found until problems with 
numerical accuracy prevent more intcrations at D=-550. 

This sequence of calculations includes a variety of behaviour and 
therefore raises a number of questions about the processes involved 
here. The first question involves the bifurcation sequence, which is 
exactly that found in the low order model of Jones et at. (1985) and 
similar to that for the third-order system examined in Tobias et at. 
(1995). It is important to discover whether the bifurcations that the 
system undergoes are simply the result of a special choice of model/ 
parameters or if the sequence is robust. The model exhibits both 
interesting temporal and spatial behaviour. Are the two linked? That 
it, does the modulation follow simply from the interaction between the 
two layers of dynamo w r aves travelling at different speeds, or is the 
spatial behaviour a secondary effect of the modulation that simply 
enhances the behaviour already present? 


6.2.2. Solutions for Constant tj 

The complicated two-layer structure for the toroidal field in the 
previous model is due to the disparity in the value of the diffusivity at 
the top and bottom of the domain. We may remove this structure by 
considering solutions to a model with constant diffusivity (i.e., a 
diffusion that is independent of both position and magnetic field 
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a) 



b) 




FIGURE 17 Chaotic solutions (a) As for Figure 13 (a), but /)= -400. (b) As or 
Figure 13(b), but D ~ — 400. (c) Chaotic attractor. The energies are now chaotically 
modulated. 


strength). For this model therefore h(z) = l and all other variables 
and parameters are as in the previous section. 

Here dynamo action to travelling waves sets in at a supercritical 
Hopf bifurcation it t £> = -83. This value is greater than that for the 
variable diflfusivity model because the diffusion is now 0(1) every- 
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FIGURE 18 Magnetic Field Spatio-temporal Chaos. As for Figure 1 1, but D = -400. 
(See Color Plate V). 

where and so the shear and a-effect must be larger to overcome 
dissipation. Again the travelling waves have constant magnetic and 
kinetic energy. Snapshots displaying the spatial dependence of the 
toroidal and poloidal field are shown in Figure 19, where it is apparent 
that the r-dependence of the solutions is less complicated than in the 
previous model. It is also interesting to note that the wavespeed of 
these solutions is c~2.8, larger than for the marginally supercritical 
solutions for the previous model so the 0(1) diffusion throughout the 
layer obviously increases the efficiency of the dynamo. This simple 
z-dependent structure is also found for the velocity perturbations. 
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FIGURE 19 Magnetic Field Travelling Waves. The figure shows the form of the 
solutions for A (top) and B (bottom) for the case where // is a constant D= -90. (See 
Color Plate VI). 


Increasing |Z)| again leads to a secondary Hopf bifurcation and 
modulated waves where the magnetic and kinetic energy both 
vacillate. The spatial structure of the modulated waves remains simple 
as shown in Figure 20, with no hint of the complicated layer structure 
and interaction between waves travelling at different speeds that 
characterised the solutions for the model with variable diffusivity. The 
modulation is clearly visible as the relative strengths of the positive 
and negative fields oscillate as the wave propagates. This is 
particularly clear in the snapshots of the poloidal field - in the first 
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FIGURE 20 Magnetic Field Modulated Waves. {// constant) The figure shows the 
form of the solutions for A (top) and B (bottom) as time progresses for D = -400. (See 
Color Plate VII). 

and fifth snapshot negative field dominates, whilst in the third the 
positive field is stronger. 

Further increase in \D\ leads to subsequent tertiary Hopf bifurcation 
to motion on a two-tours in 'energy-space' (shown in Fig. 21) which 
breaks down to form the chaotic attractor. For these solutions the 
spatial structure is complicated. The nodes are now irregular and 
breathe erratically, as well as travelling quickly to the right. The 
velocity field behaves in much the same way wild oscillations are 
superimposed on a net propagation of the waves to the right. If the 
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FIGURE 21 Two-torus and chaotic attractor (// constant): Plots of the trajectories in 
the ( Eh , FA) plane for (a) D 1000. (b) D — 1500. 


dynamo number is increased further these oscillations become larger in 
amplitude. 

This model therefore undergoes the same sequence of bifurcations 
as the more realistic model discussed earlier. The travelling waves 
become unstable in a series of Hopi bilurcations leading to motion on 
a torus which breaks down, leaving a chaotic attractor. The 
modulation is therefore a direct result of including the Malkus- 
Proctor mechanism and the equation for the large-scale velocity 
perturbations. Any complicated spatial behaviour is a secondary 
effect, which is only present because of the bifurcations to modulated 
waves. (If travelling waves remain stable then the spatial structure 
stays simple, as in the local quenching models.) Whilst the spatial 
dependence may enhance the temporal behaviour and encourage the 
transition to spatio-temporal chaos, it is not an integral part ol the 
mechanism. 
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6.2.3. Changing the Magnetic Prandtl Number 

We now return to the initial Parker model with diffusivity varying 
across the layer, as this model is more relevant to the generation of 
fields at the base of the convection zone. The consequence of varying 
the magnetic Prandtl number on the nature of solutions is investigated. 
In the low-order model studied by Weiss et al. (1984) and Jones et at. 
( 1 985), the bifurcation sequence was fairly sensitive to the choice of the 
value of the magnetic Prandtl number with the periodic solution 
remaining stable unless r< 1. Clearly if the same is true for the PDEs 
investigated here, then this would have serious consequences for the 
claim that the bifurcation sequence and route to chaos found here are 
robust and therefore astrophysically relevant. 

For this reason, results are now reported for r = 0.5 and r = 2. 
(Changing the magnetic Prandtl number by an order of magnitude 
presents problems for the numerical scheme though this is addressed in 
the next section.) The results for r = 0.5 are similar to those for r = 1. 
The bifurcation sequence is the same and the only difference is in the 
value of the dynamo number at which the bifurcations take place. The 
bifurcations for this model occur earlier (at lower values of |£>|) than 
lor when r = 1 . This is to be expected as the smaller the value of r, the 
larger the effect of a given magnetic field strength on the velocity 
perturbations. The time series for the magnetic energy for three types 
of solution (travelling waves, modulated waves and quasiperiodic 
waves) are shown in Figure 22. 

Not only is the bifurcation structure the same for this parameter 
value, but the spatial dependence of solutions is similar to that for 
r = 1 - the two-layer structure of the toroidal field remains, and the 
four nodes of magnetic field persist in the lower layer. The spatial 
dependence is therefore the same as that shown in Figure 14. 

For r = 2 a similar bifurcation sequence is again observed. This time 
the bifurcations are delayed in comparison with the r = 1 case. Here 
though the spatial dependence of the solutions is slightly different. The 
series of snapshots for the toroidal field for values of the dynamo 
number that give modulated waves in Figure 23, shows that although 
there are still two layers of magnetic field, the slow-moving magnetic 
field (at the bottom of the box) has a similar horizontal scale to that in 
the fast-moving magnetic field. This can be ascribed to the increase in 
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c) 



FIGURE 22 Solutions for r = 0.5. Magnetic energies for (a) D = —50, (b) D ™ -200, (c) 
£> = -500. 

r leading to a larger scale for the magnetically driven velocity per- 
turbations. Apart from the scale of the field and velocity perturba- 
tions, the solutions are similar to those for r= 1, and the sequence of 
bifurcations is certainly the same - a succession of Hopf bifurcations 
followed by breakdown of a tours leading to spatio-temporal of 
chaos. The spatial dependence of the chaotic solutions is slightly 
different to that for r=l, but in all other respects the solutions are 
similar. 


6,2 A. A Note on Very Low Prandtl Numbers 

So far in this section only 0(1) values of r have been considered. As 
noted in the previous section, changing the value of the magnetic 
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FIGURE: 23 Magnetic field Modulated waves for t 2: As for Figure 14, hut r 2 
and /) ~ 250. This time the magnetic field generated by the velocity perturbations has 

the same horizontal scale as that generated by the imposed shear. (See Color Plate VIII). 


Prundtl number by an order of magnitude presents problems for the 
numerical scheme in the model where the diffusivity is a function of 
position. However the effects of very low Prandtl numbers can be 
examined in the constant-// model discussed in Section 6.2.2. A 
thorough investigation of parameter space is certainly not possible in 
this case, as it takes very many magnetic diffusion times for transients 
to die away. For this reason only some preliminary results have been 
obtained. The results of this section have led to investigation of the 
consequences of including the Malkus-Proctor effect with a low value 
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of r in a global dynamo model where the variation of parameters with 
latitude is included (Tobias 1996a, 1997). 

Here the magnetic Prandtl number is set to be 0.01 and numerical 
integrations are performed for a selection of dynamo numbers. We are 
primarily interested in the behaviour of solutions for large(ish) values 
of \D\. The models already discussed do produce chaotic modulation 
of the basic magnetic cycle, but this is a weak effect, with the 
amplitude of modulation small in comparison to the amplitude of the 
basic cycle. So it is important to ascertain whether the low-r limit 
yields solutions where the magnetic energy is substantially modulated 
and where periods of greatly reduced activity (corresponding to 
minima) can be detected. 

The first value of the dynamo number investigated was £> = -2100, 
which is far in the nonlinear regime. The results are shown in Figure 24. 
The times series for both the magnetic and kinetic energy are periodic. 
This indicates that the solutions have undergone both a primary Hopf 
bifurcation to travelling wave (constant energy) solutions and a 
secondary Hopf bifurcation to quasiperiodic solutions (with periodic 



FIGURE 24 Period-2 orbit: Low r solutions for r = 0.0L /)=-2100. (a) Average 
magnetic energy as a function of time. It is important to remember that the basic 
magnetic cycle has been removed from this time series by averaging. The time series 
therfore gives the amplitude of the modulation, there are episodes when the magnetic 
energy reduces drastically to zero, (b) As for (a), but for the average kinetic energy, (c) 
Phase-plane trajectory inergy space. Graph of ( B ~) vs (v“). 
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average energies). Apparently the solution has also undergone a 
subsequent period-doubling bifurcation. The figure clearly shows that 
the modulation of the basic cycle (which has been removed by the 
averaging procedure) is a large effect and that there are episodes when 
the magnetic energy is very small. 

Chaotic modulation can be found if the magnitude of the dynamo 
number is increased further, as in Figure 25 where solutions for 
D = - 2500 are shown. The kinetic and magnetic energy are both 
irregularly modulated and the attractor in phase space appears to be 
chaotic. The modulation is still a large effect and many episodes of 
reduced activity are detected here. 

These results give an indication of how solutions behave if the 
magnetic Prandtl number is reduced by an order of magnitude. It 
seems a though the low-r limit leads to significant amplitude 
modulation of the basic magnetic activity. The low value of the 
Prandtl number makes larger velocity perturbations easier to excite, 
thereby increasing the amplitude of the modulation. 



FIGURE 25 Chaotic orbit: Lot r solutions for r = 0.01, £>--2500. (a) Average 
magnetic energy as a function of time, (b) As for (a), but for the average kinetic energy, 
(c) Phase-plane trajectory in energy space. Graph of { b 2 ) vs (v 2 ). Both the magnetic and 
kinetic energy are now chaotically modulated. 
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6.3. Summary for Malkus-Proctor Nonlinearities 

The coupling of the dynamo equations to the mean momentum 
equation has been studied in two separate cases. In the first problem, 
velocity perturbations that were independent of z were driven by the 
averaged Lorentz force. This system produced travelling wave 
solutions that were stable for all values of the dynamo number D 
and in this respect the solutions were the similar to those found for the 
simple quenching mechanisms discussed in Section 5. It is interesting 
to note that the magnetic field generated velocity perturbations away 
from the region of strongest magnetic field, differing from the 
quenching mechanisms where only local action of the magnetic field 
can be modelled. The wavespeed increased substantially with 
increasing dynamo number. 

Interesting time-dependent behaviour was only achieved if the local 
field was allowed to interact with the velocity in the mean momentum 
equation. The important features of the bifurcation sequence are 
robust to changes in the parameters and even to changes in the profile 
for the diffusivity. The trivial state becomes unstable to dynamo action 
in a Hopf bifurcation (which is usually supercritical) to travelling wave 
solutions. Solutions are steady in a frame moving with the wavespeed 
and the average magnetic and kinetic energy are constant. These 
travelling waves lose stability in a secondary Hopf bifurcation as the 
dynamo number is further increased - the average energies now ; 
oscillate about a mean value. A further (tertiary) Hopf bifurcation 
leads to motion on a two-tours in "average energy space' (or a three- 
torus in ‘magnetic field space 1 ). Solution that exhibit spatio-temporal 
chaos are found as the dynamo number is increased further and the 
torus breaks down. It is also possible to find other bifurcations 
(pitchfork and gluing bifurcations) if the parameters are chosen 
carefully, but these are not essential parts of the bifurcations to chaotic 
solutions. Complicated spatial dependence of solutions is also found 
but this is a secondary effect. The most significant modulation was 
obtained when the magnetic Prandtl number (r) was set to be small. 

The route to chaos is therefore exactly the same as that found in the 
sixth-order model of Weiss et al. (1984). In coordinates moving w'ith 
the wavespeed, the initial Hopf bifurcation becomes a pitchfork and 
the bifurcation structure is now the same as the simpler model of Jones 
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ct ciL (1985). It is important to understand why the system must 
undergo three Hopf bifurcations before chaotic behaviour can be 
found. The frequency introduced in the Hopf bifurcation to travelling 
waves is removable (as described above) and if this frequency is 
removed the travelling waves can be thought of as steady solutions. 
This was noted by Ruelle (1973) for waves in a rotating system (but the 
extension to periodic systems is immediate). With the frequency 
removed, the fixed point must undergo two further Hopf bifurcations 
before frequency-locking and the formation of a chaotic attractor via 
the breakdown of a torus can occur (Rand, 1982). So here the route to 
chaos is slightly different from that for the case where the system has 
boundaries in Tobias (1996a). There the initial Hopf bifurcation was 
not to travelling waves, but to genuinely periodic solutions and the 
frequency of oscillation was not removable. Therefore only one more 
Hopf bifurcation was required to introduce the possibility of 
frequency-locking and breakdown of the torus. 

It is therefore only due to the artificial nature of the periodic 
boundary conditions that the system must undergo a tertiary 
bifurcation before chaotic solutions may be found. In systems with 
more realistic boundary conditions this may occur after only two Hopf 
bifurcations and this is therefore the genetic scenario for transition 
to chaos in these dynamo models. 


7. CONCLUSION 

The explanation of the origin of magnetic activity in the sun and stars 
provides a challenge to the theoretician. In the past this challenge has 
been met in a rather ad-hoc and haphazard way. Many nonlinear 
dynamo models have been proposed and results exhibited, but relatively 
little has been learned of the underlying structure of nonlinear dynamos. 
The models have only scraped the surface of the problem - the results 
are clear, but what are their consequences for stellar magnetic field 
generation? Is the behaviour found model-dependent, and if not why is it 
a robust feature of the dynamo problem? 

By carrying out a systematic investigation of a simple dynamo 
model, I have attempted to provide a clear and coherent picture of the 
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role of nonlinear effects in mean-field dynamos. Given that there is so 
much that is poorly understood about the nature of the dynamo 
problem, it is only by interpreting the results in the light of the recent 
advances in dynamical systems theory that progress may be made. We 
need to separate out the robust behaviour and discard any results that 
are model-dependent. An understanding of the fundamental proper- 
ties of nonlinear systems is essential in order to achieve this. 

Local travelling solutions do yield a greater understanding of the 
consequences of nonlinearities for the radial dependence of solutions. 
Moreover, much can also be learned about the consequences of 
including simple quenching mechanisms. Quenching of the turbulent 
diffusion is found to have a more subtle effect than either a or w- 
quenching, as it modifies the spatial structure of solutions. However, 
the unrealistic choice of ( periodic) lateral boundary conditions lor this 
local model implies that any global properties (e.g., wavespeed, 
frequency) of such a dynamo may be misleading (Tobias ct ai . , 1997). 
This is an important result, as it demonstrates that over-simplified 
models may lead to incorrect conclusions about the behaviour of 
dynamos. 

Quenching mechanisms, although simple to implement, are un- 
realistic parameterizations of the back-reaction of the magnetic field 
on the velocity. These mechanisms are instantaneous in time and local 
in space a severe restriction on the possible consequences of 
including nonlinearities. It is interesting to note that, because of the 
simplicity (and local nature) of the quenching mechanisms often used 
in nonlinear mean field dynamos, the behaviour of the solutions is not 
critically dependent on the precise form of the quenching chosen. This 
was shown to be the casse as the imposition of both ‘averaged’ (one- 
dimensional) and ‘local’ (two-dimensional) quenching mechanisms 
yielded similar results. 

Increasingly complicated spatio-temporal behaviour may be 
achieved in a natural way by including an equation for the mean 
velocity perturbation driven by the magnetic field (the Malkus-Proctor 
effect) as in Section 6. This nonlinearity is more complicated than the 
simple prescribed quenching formulae for two reasons. The first is that 
this nonlinearity is non-local in space. The velocity perturbations are 
not driven where the field is largest but where the gradient of the 
poloidal and toroidal fields interact. Hence, for this nonlinearity, the 
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results are sensitive to the choice of averaged or local interaction. 
Averaged interaction yielded results similar to those for the quenching 
mechanisms already discussed, with travelling waves remaining stable. 
However, if local interactions are allowed then the system displays 
very different behaviour with the travelling waves becoming unstable 
to solutions with more complicated time-dependence. 

Moreover the torsional oscillations observed in the sun can be 
attributed to such a velocity perturbation. The addition of a separate 
evolution equation then allows for a time-lag between the generation 
of the field and its effect on the velocity field. This time-lag is 
controlled by the magnetic Prandtl number r. It is important to note 
that the route to chaos is via a succession of Hopf bifurcations and the 
breakdown of a torus — similar to that for the low-order model of 
Tobias et al. (1995). This is a very important result as it implies that 
this mechanism is a robust feature of including the interaction of the 
magnetic field with the velocity in a realistic way. 

What do the above results imply for modelling dynamo action in the 
sun and other stars? By studying the role of nonlinearities in isolation, 
one can interpret the results of numerical simulations which include 
latitudinal dependence of parameters and realistic boundary condi- 
tions. Indeed the results included in this paper have already been used 
(Tobias 1996a, 1997) to construct global models of the solar dynamo. 
In those models the 1ow-t limit was successfully exploited to 
investigate the nature of Grand Minima, and to ascertain the relative 
importance of modulational mechanisms (i.e., the Malkus-Proctor 
mechanism and Parity interactions) on solar dynamo models. The 
results are promising but further investigations are required to 
understand the nature of the interaction of turbulent convection and 
magnetic fields so that more realistic models may be constructed in the 
future. 
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